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A comparative numerical study of meshing 
functionals for variational mesh adaptation* 




Weizhang 


Abstract 


We present a comparative numerical study for three functionals used for variational mesh 
adaptation. One of them is a generalization of Winslow’s variable diffusion functional while 
the others are based on equidistribution and alignment. These functionals are known to 
have nice theoretical properties and work well for most mesh adaptation problems either as 
a stand-alone variational method or combined within the moving mesh framework. Their 
performance is investigated numerically in terms of equidistribution and alignment mesh 
quality measures. Numerical results in 2D and 3D are presented. 
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1 Introduction 

Variational mesh adaptation is an important type of mesh adaptation method and has received 
considerable attention from scientists and engineers; e.g., see the books da ESI ESI El and 
references therein. It also serves as the base of a number of commonly used adaptive moving 
mesh methods (e.g., see [3113111 El). In the variational approach, an adaptive mesh is 
generated as the image of a reference mesh under a coordinate transformation and such a 
coordinate transformation is determined as a minimizer of a certain meshing functional. A 
number of meshing fnnctionals have been developed in the past (cf. the above mentioned books). 
For example, Winslow |25) proposed an equipotential method based on variable diffusion. 
Brackbill and Saltzman [ 3 ] developed a method by combining mesh concentration, smoothness, 
and orthogonality. Dvinsky |B] used the energy of harmonic mappings as his meshing functional, 
while Knupp [21 and Knupp and Robidoux m developed functionals based on the idea of 
conditioning the Jacobian matrix of the coordinate transformation. More recently, Huang [ 7 ] 
and Huang and Russell m proposed functionals based on the so-called equidistribution and 
alignment conditions. 
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With variational mesh adaptation, the mesh concentration is typically controlled through a 
scalar or a matrix-valued function, often referred to as the metric tensor or monitor function 
and defined based on some error estimates and/or physical considerations. While most of 
the meshing functionals have been developed with physical or geometric intuitions and have 
various levels of success in the adaptive numerical solution of partial differential equations 
(PDEs) and other applications, there is only a limited understanding on how the metric tensor 
affects the behavior of the mesh. An attempt to alleviate this lack of understanding was 
made by Cao et al. |1] for a generalization of Winslow’s variable diffusion functional. They 
showed that a significant change in an eigenvalue of the metric tensor along the corresponding 
eigendirection (first increasing and then decreasing, or vice versa) will result in adaptation 
of coordinate lines along this direction, although this adaptation competes with far more 
complicated effects, including those from changes in eigenvectors and other eigenvalues and 
the effects of the shapes of the physical and computational domains and the mesh point 
distribution on the boundaries. In 0115], two functionals have been developed based directly 
on the equidistribution and alignment conditions. These two conditions provide a complete 
characterization of the mesh elements through the metric tensor [ 7 ] . Minimizing the functionals 
leads to meshes which tend to satisfy the conditions in an integral sense. Nevertheless, this 
characterization is only qualitative, and how closely the mesh satisfies the conditions depends on 
the boundary correspondence between the computational and physical domains and the mesh 
point distribution on the boundaries. Thus, numerical studies, especially comparative ones, 
are useful, and often necessary, in understanding how the mesh adaptation for those meshing 
functionals is controlled precisely by the metric tensor. There do exist a few comparative 
numerical studies for meshing functionals. For example, a gallery of (adaptive and non-adaptive) 
meshes is given in m for a number of meshing functionals. Some comparative meshes are 
given in m for the harmonic mapping functional [Bj and the subsequent functional based on 
equidistribution and alignment [ 7 ] . 

The main objective of this work is to present a comparative study for three of the most 
appealing meshing functionals, a generalization of Winslow’s variable diffusion functional 
(cf. Q) and two functionals based on equidistribution and alignment (cf. ( |11[ ) and ( |13[ )). They 
are selected because (§ and (0 have been known to work well for many problems (e.g., 
see [1121 El US] [HISS]) while ( |I^ is similar to ( [IT] ) and has some very nice theoretical properties 
(cf. §3.2). Another motivation is to present some three dimensional numerical results for 


those functionals. Critical for our study is to perform the substantial computations using the 
improved implementation of the variational methods introduced in m- In a sharp contrast to 
the situation in two dimensions, very little work has been done with variational mesh adaptation 
and adaptive moving mesh methods in three dimensions (e.g., see [151 E5|). It is particularly 
interesting to see how the functionals perform in this case. 

An outline of the paper is given as follows. We describe the basic ideas of the variational mesh 


adaptation and its direct discretization (that is, first to discretize and then optimize) in §2 In §3 


we introduce the three functionals to be studied for the numerical comparison, a generalization 
of Winslow’s variable diffusion functional and two functionals based on equidistribution and 
alignment. Numerical results and example adaptive meshes are given in [§4] followed by 
conclusions in §5 
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2 Variational mesh adaptation 


In variational mesh adaptation, an adaptive mesh is generated as the image of a reference mesh 
under a coordinate transformation. Denote the physical domain by D C {d> 1), and assume 
that we are given a computational domain Dc C and a quasi-uniform mesh The thereon (in 
this work we consider only simplicial meshes). In many situations we can choose Dc to be the 
unit square/cube or simply D. Denote the coordinate transformation by a; = x {^): Dc —>■ D 
and its inverse by ^ = ${x): D —?■ Dc. Such a coordinate transformation (more precisely, its 
inverse) is determined as a minimizer of a meshing functional. Most of the existing meshing 
functionals can be cast in a general form as 

I[^]= f G(J,det(J),M,a;)da;, (1) 

Jn 

where G is a smooth function, J is the Jacobian matrix of ^ = ^{x), det(JI) denotes the 
determinant of J, and M = M(a;) is the metric tensor supplied by the user to control mesh 
concentration. We assume that M = M(a;) is a symmetric and uniformly positive definite d-hy-d 
matrix-valued function on D. Notice that (Q is formulated in terms of the inverse coordinate 
transformation. One reason for this is that this form is less likely to produce singular meshes [B] . 
Another reason is that M is a function of x and thus finding the functional derivative of /[^] 
will not directly involve the derivatives of M. This is convenient since in practice M is known 
only at the vertices of a mesh and its derivatives are not cheap to find. The main disadvantage 
of the formulation in this form is that ^ = ${x) (or its numerical approximation) does not give 
the physical mesh directly. This is remedied either by interchanging the roles of the independent 
and dependent variables in the Euler-Lagrange equation of I[^] (e.g., see [15]) or, in a recently 
developed implementation (see below), computing the new physical mesh from a computational 
one using linear interpolation. 

A minimizer of ([^ can be found numerically in the MMPDE (moving mesh PDE) framework. 
A conventional implementation m is to find the functional derivative of Q and then define 
the MMPDE as the gradient flow equation of the functional. Having been transformed by 
interchanging the roles of the dependent and independent variables, the MMPDE can be 
discretized on The and a system of equations for the nodal velocities is obtained. Einally, the 
new mesh is obtained by integrating the mesh equation over a time step. 

A much simpler implementation was proposed recently by Huang and Kamenski m- Instead 
of utilizing the MMPDE directly, the new implementation first discretizes the functional on the 
current mesh 7h and then, following the idea of the MMPDE, dehnes the mesh equation as the 
gradient equation of the discretized functional (with respect to the computational coordinates 
of the vertices). To be specific, denote by Xi, and z = 1,..., W the coordinates of the 
vertices of the current physical mesh (Th), the reference mesh {Thc^, and the computational 
mesh {Thc)j respectively. We assume that these meshes have the same numbers of the elements 
and vertices and the same connectivity. For any element K £ Th (with vertices xf, z = 0,..., d), 
the corresponding element in The is denoted by Ke (with vertices z = 0,..., d). The edge 
matrices for K and Ke are defined as 

Ek = [x^ -x^,...,x^ -x^], Ek, = [^1 -^o]- 


Let uji be the element patch associated with vertex Xi (i.e., the collection of the elements 
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containing Xi as a vertex). Then, the equation for the nodal velocities reads as 


,K 


i - 1 , . . . , Ny , 


(2) 


where |iir| is the volume of K, is the local mesh velocity associated with vertex Xi in K, 
ix denotes the local index of ajj in iT, r > 0 is a constant parameter used to adjust the time 
scale of mesh movement, and P = (Pi,..., Pnv) is a positive function used to make the mesh 
equation to have desired invariance properties. Although the parameter r can be absorbed in 
P, using two parameters has the advantage that the role of each parameter is clear: r for the 
time scale of mesh movement while P for the invariance properties. Ideally, r should be chosen 
such that the mesh movement has the same scale as the physical equation. Unfortunately, there 
is no theoretical analysis for this yet and trial and error is still the most practical way to choose 
r. Numerical experience shows that a value in the range [0.01, 0.1] seems to work well for most 
problems m- In our computation, we choose P such that the equation is invariant under 
the scaling transformation M —?■ cM for all non-zero constants c (cf. §3): in variational mesh 
adaptation it is the relative distribution of M over the physical domain (instead of its absolute 
distribution) that determines the variation of the mesh density and therefore it is essential for 
the moving mesh equation to be invariant under the scaling transformation of M. 

The local velocities are given by 


(uf)^ 




= -E 


idG 


dG 


K 


det {EkJ 1 


9det(JI) det(Px) 




K 


= -E- 


where the derivatives of G with respect to JJ and det(JI) 
evaluated as 


(see 0, @, and @ below) 


(3) 


are 


dG dG / dePExJ ^ 

AT = AT \^k,Ek\ 


dG 


dG 


Ek,E, 


det (Ek) 

i det(Pii-J 


,M{xk),xk 


9det(J[) 9det(JI) v = ^ ’ det(Pii') 

The second equation in 0 is an inherent property resulting directly from the differentiation 
of the discretized meshing functional; it states that the centroid of K stays fixed if only the 
contribution from K is taken into account. 

The above mesh equation should be modified properly for boundary vertices. For example, if 
is a fixed boundary vertex, we replace the corresponding equation by 

t-- 

When is allowed to move on a boundary curve/surface represented by 


m = 0, 

then the mesh velocity ^ needs to be modified such that its normal component along the 
curve or surface is zero, i.e., 

= 0 . 
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Mesh equation Q defines the movement of the nodes of the computational mesh The starting 
from the reference mesh The ^n- The equation can be integrated in time to obtain the 
computational mesh at tn+i- For notational simplicity, we denote the computational mesh at 
tn+i by The well. Notice that the physical mesh Th is fixed during the time integration from 
tn to tn+i- Meshes The and Th define a correspondence 

Th = ^{The)- 


The new physical mesh is computed by means of linear interpolation as 

fh = ^{fhe), 


where The is the reference mesh on ^}e■ The computational mesh plays the role of an intermediate 
variable. 

Recall that the mesh concentration in variational mesh adaptation is controlled through the 
metric tensor M = M(a:). Such a metric tensor can be defined based on physical or geometric 
considerations or some error estimates. For example, for the norm of the error of piecewise 
linear interpolation on simplicial meshes, the optimal metric tensor [HKTB] (also see m (5.192)]) 
is 

M = det (a/ + \H{u)\)~^ [al + \H{u)\] , (4) 

where H(u) is the Hessian of function u, \H{u)\ is the eigen-decomposition of H{u) with the 
eigenvalues being replaced by their absolute values, and the regularization parameter a > 0 is 
chosen such that 


det(M)2(i®= / det {al + \H{u)\)‘‘+* dx = 2 / det dx. 


jn JQ jn 

In practical computation, u is typically unknown, and only the approximations to its values at 
the vertices are available. For this reason (and even in situations where an analytical expression 
for u is available), the Hessian in (Q is replaced by an approximation obtained by a Hessian 
recovery technique from the nodal values of u or the approximations of the nodal values of u. 
A number of such techniques are known to produce nonconvergent recovered Hessians from 
a linear finite element approximation (e.g., see Kamenski [E])- Nevertheless, it is shown by 
Kamenski and Huang m that a linear finite element solution of an elliptic BVP converges 
at a second order rate as the mesh is refined if the recovered Hessian used to generate the 
adaptive mesh satisfies a closeness assumption. Numerical experiment shows that this closeness 
assumption is satisfied by the approximate Hessian obtained with commonly used Hessian 
recovery methods. We use a Hessian recovery method based on a least squares fit: a quadratic 
polynomial is constructed locally for each vertex via least squares fitting to neighboring nodal 
function values and an approximate Hessian at the vertex is then obtained by differentiating 
the polynomial. 


3 Meshing functionals 

Here we introduce the three meshing functionals used in the numerical study. A generalization 
of Winslow’s variable diffusion functional and the two functionals based on equidistribution 
and alignment are selected because they are reasonably simple, have nice theoretical properties, 
and are known to work well for many problems. 
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3.1 Winslow's functional based on variable diffusion 


The first functional is the variable diffusion proposed by Winslow [25] . It uses the system of 
elliptic PDEs 

-V • = 0, i = 

for generating adaptive meshes, where w = w{x) > 0 is the weight function. This system 
mimics a (steady-state) diffusion process with a heterogeneous diffusion coefficient w{x). It is 
the Euler-Lagrange equation of the functional 

= d f = d f w{x)tT{M'^)dx, (5) 

^ JO ^ J Jo 

where tr(-) is the trace of a matrix. A generalization of this functional to allow a diffusion 
tensor is 

= X / tr(JIM“^Jl'^) dx. (6) 

J Jo 

This functional has been used by a number of researchers, e.g., see Huang and Russell caiii], 
Li et al. [22] j and Beckett et al. |2] . It is coercive and convex m Example 6.2.1]. Thus, under a 
suitable boundary condition (such as the Dirichlet boundary condition with d^c being mapped 
onto (9H), the functional (§ has a unique minimizer. 

For this functional, the derivatives of G with respect to J and det(JJ) needed in Q are 

{G= itr(JM-ij^), 

f = M-^F, (7) 

i Sdt?(J) = O' 

The interested reader is referred to m for the derivation. 

The balancing function in ([^ is chosen as P = det(M)d. With this choice, ([^ is invariant 
under the scaling transformation M —)■ cM. 


3.2 Functionals based on equidistribution and alignment 


The other functionals are based on the equidistribution and alignment conditions. These 
conditions provide a full mathematical characterization of a non-uniform mesh. Indeed, any 
non-uniform mesh can be viewed as a uniform one in the metric specified by a tensor. Moreover, 
a mesh is uniform in the metric specified by the metric tensor M = M(a;) if and only if it satisfies 
the equidistribution and alignment conditions associated with M [iniiisi. In the continuous 
form, they are 


where 


equidistribution: 

alignment: 


det(JI) ^det(M )2 = Va; e 0 

\^^c\ 

^ tr(JM-^= det(JM-^s, V® G H, 


cr = [ det(M) 2 (ia;. 
Jo 


( 8 ) 

(9) 

( 10 ) 
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These conditions require the mesh elements to have the same size (equidistribution) and be 
equilateral (alignment) in the metric M, respectively. The alignment condition also implies that 
the elements are aligned with M in the sense that the principal directions of the circumscribed 
ellipsoid of each element coincide with the eigen-directions of M while the lengths of the principal 
axes of the ellipsoid are reciprocally proportional to the square roots of the eigenvalues of M. 
The first functional based on equidistribution and alignment, proposed in 13 , is 

I[$]=eJ ^det(M) (tr(JM-ij^)) dx + {1 - 2e)d^ J ^d^t^ dx, (11) 

where 9 G (0,1) and p > 0 are dimensionless parameters. Loosely speaking, the first and second 
terms correspond to the equidistribution and alignment conditions, respectively. The terms are 
dimensionally homogeneous and the balance between them is controlled by the dimensionless 
parameter 9. For 0 < 9 < ^, dp > 2, and p> 1, the functional is coercive and polyconvex and 
has a minimizer m Example 6.2.2]. Moreover, for 9 = ^ and dp = 2 it reduces to 

[ \/det(M) tr(JIM“^Jl'^) dx, 

which is exactly the energy functional of a harmonic mapping from to Qc (cf- IS])- 
For the functional , we have 

G = gv/det(M)(tr(JM-^F))^ + (1 - 2g)d^ v/det(M) , 

' ^ = dp0v'det(M)(tr(JM-iF))^“^M-iF, 

a^j) =P(1 -26l)d^ det (M)^det (JI)^“^ 

In the computation, we use {p, 0) = (2, ^). p = 2 is the smallest integer to satisfy dp >2 for 
d = 1, 2, and 3. The choice of 0 = 1/3 is in the range (0,1/2] for the functional to be polyconvex 
while giving a bigger weight to the equdistribution condition. This set of the values works well 
for all tested problems. The balancing function in (|^ is chosen to be P = det(M)^^, so that 
([^ is invariant under the scaling transformation M —)■ c M. 

The second functional based on equidistribution and alignment is 


m = 01 / Jde 




dx 


+ 


dp(d-2) r Ir-, dp ^ dp / ^ \ 

02d2(d-i) / det(M)5^^“J^Met(J)^ (tr(J-^MJ-i)j d® 


+ 


9a f 

o-P+^ Jn 


^ ^ V Vdet(M) _ 


dx, (13) 


where p > 1, > 0, and 0* > 0 (i = 1,..., 4) are parameters. The first three terms in (13) are 

dimensionally homogeneous in M and J while the last term has the same dimension in M as the 
other terms. This functional was proposed in [T51 (6.120)] to avoid singularity of the coordinate 
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transformation. Indeed, if 9^ — 6 i — O 2 > 0, then the functional is coercive and polyconvex and 
has a minimizer satisfying det(JI) > 0 in jl51 Example 6.2.3]. 

In the computation, we choose 61 = 62 = 5 , 03 = 1, 04 = O.l,p = 2, v = 1, and the balancing 




function P = det(M)~. These choices are based on the functional (11) and the desire to keep 
the fourth term relatively small. 

For this functional, we have 


G = ^i^de 


tr 


-IF 


dp 

~2 


90 

aji 


ap\a—A) \ dp \ dp / 

+02c( det(M)2i d-il (let(JI) <^-1 ^tr(JI 

+(03 - 01 - 02)d'^ (v/det(M)) 

+ (^/d<M))'^"det(J)-^ 

^ ^ ^_ 1 

= 0 idp\/det(M)(tr(JM-iF)) " M-^F 


det(J)P 


(14) 


dp 

-In \ 2(d-l) 


-^d 2V1) det(M)5(i-^) det(J)^ (tr(JJ-^ 

det(M)5(^-^Fet(J)^-^(tr(J-^MJ-i)) 
+(03 — 01 — 92)pd^ (^-ydet(M)^ ^det(JI)P-^ 


dp 

2{d-l) 


641/ 

crP+^ 


^A/det(M)^ det(JJ) 


-u-l 


4 Numerical experiments 

In the following we consider a number of examples in two and three dimensions. For a given 
function we consider M defined in (|^ which is optimal for minimizing the norm of the linear 
interpolation error of this function and compare meshes obtained from using the functionals of 
Winslow ([^ (W), Huang ( [IT| ) (H), and Huang and Russell (13) (HR). 

To assess the quality of the generated meshes, we compare the norm of the linear 
interpolation error and the equidistribution and alignment mesh quality measures, which 
describe how far the mesh is from being uniform in the metric defined by M. The element-wise 
quality measures are based on ^ and and defined as 


Qeq,K — 


det(JI/i') ^det(Mic)2 


Qali,K — 


tr () 

cidet(JIii-M]^^JI^)^ 


(15) 


while for the overall mesh quality measures we take their root-mean-squared values. 


Qeq — 




1 

N 


E 


K&Th 


eq^K ’ 


Qali — 




1 

N 


E Q‘ili,K- 


(16) 


K&Th 


If the mesh is uniform with respect to M, then Q^q = QaH = 1; if the mesh is far from being 
uniform with respect to M, then Qeq and Qaii will become large. In other words, these quality 
measures describe how well the volume (measured by Qeq) and the shape and orientation 
(measured by Qaii) of mesh elements correspond to the desired size and shape prescribed by M 
(see [5] for more details on the mesh quality measures). 
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4.1 Two dimensions 


First, we consider two dimensional meshes constructed for the unit square 17 = (0,1) x (0,1) 
and the test functions 


Example 4.1. 


u = tanh (—30 [y — 0.5 — 0.25 sin(27rrE)]) 


Example 4.2. 

u = tanh(25y) — tanh(25(x — y — 0.5)). 

Example meshes, close-ups, as well as the mesh quality measures and the interpolation 
error are given in Figs. [^and[^ 

For these examples, all three functionals provide good size and shape adaptation. A closer 
look at the mesh quality measures shows that, although all three functionals provide comparable 
meshes which are reasonably close to the prescribed metric tensor, meshes constructed using 
H and HR functionals have better correspondence to the prescribed metric tensor. In both 
two-dimensional examples, H and HR functionals provide very similar grids which are closer to 
the prescribed size and shape (that is, smaller values of Qeq and Qaii)- This is also reflected in 


the error of the linear interpolation: HR functional (13) provides the smallest error, followed by 
H functional ( [IT| ) and then W functional ([^. It seems that W functional is a bit too aggressive 
in moving nodes towards the neighborhood of areas of interest, providing a higher density of the 
nodes along the anisotropic features of the given function while coarsening out the mesh nearby, 
leading to a steeper element size gradation. Interestingly, for both examples the convergence of 
the linear interpolation error for W functional (Figs. and slows down near N = 10^ and 
returns to the order 0{N~^) as the mesh is refined {N is the number of mesh elements). It is 
unclear to us what causes this for W functional. 

For Example |4.1| we also compute adaptive meshes for the metric tensor which is optimal 
for the semi norm error (see, e.g., m for details on the metric tensor). The results shown 
in Fig. [^are essentially the same as in Fig. HR functional (13) provides the smallest error, 
followed by H functional (11) and then W functional For this metric tensor, the error for 
the H and HR functionals (Fig. is slightly larger than in Fig. which is not surprising, since 
the metric tensor chosen for Fig.[^is optimal for the error. Thus, adding the equidistribution 
property to the meshing functional seems to help to obtain meshes which are closer to fulfilling 
the prescribed properties. Interestingly, the error for the W functional seems to be a bit 
smaller if the Ef^ metric tensor is used. This can be explained by the fact that the W functional 
does not incorporate the equidistribution property and, thus, doesn’t exactly generate a mesh 
which is uniform with respect to the provided metric tensor. Thus, it is not quite clear if one is 
able to generate the optimal mesh when using the W functional. 


4.2 Three dimensions 

In three dimensions, we consider the unit cube H = (0,1) x (0,1) x (0,1) and the following test 
functions. 
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(d) Qeq VS. N 



(a) Winslow’s ^ 




(c) Huang and Russell’s (131 
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(e) Qaii vs. N 




Figure 1: Example 


4.1 


example meshes (left), close-ups near the wave tip (middle) and in 


the middle (right), mesh quality measures, and interpolation error (black line 
represents A^“^). 
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Figure 2: Example 4.2 example meshes (left), close-ups near the wave meeting the boundary 


layer (middle) and in the right bottom corner (right), mesh quality measures, and L? 
interpolation error (black line represents N~^). 
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(a) Winslow’s 
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(b) Huang’s (111 







Figure 3: Example 4.1 using a metric tensor for the semi-norm: 


close-ups near the wave tip (middle) and in the middle (right) 
and interpolation error (black line represents A^“^). 


example meshes (left), 
mesh quality measures, 
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Example 4.3. 


u = tanh (^30 (4x — 2.0)^ + {Ay 
+ tanh (^30 {Ax — 2.5)^ + {Ay 
+ tanh ^30 {Ax — 2.5)^ + {Ay 
+ tanh (^30 {Ax — 1.5)^ + {Ay 
+ tanh (^30 {Ax — 1.5)^ + {Ay 
+ tanh ^30 {Ax — 2.5)^ + {Ay 
+ tanh (^30 {Ax — 2.5)^ + {Ay 
+ tanh (^30 {Ax — 1.5)^ + {Ay 
+ tanh ^30 {Ax — 1.5)^ + {Ay 


2.0)2 + {Az 

2.5) 2 + {Az 

1.5) 2 + {Az 

2.5) 2 + {Az 

1.5) 2 + (4z 

2.5) 2 + {Az 

1.5) 2 + {Az 

2.5) 2 + {Az 

1.5) 2 + {Az 


2 . 0)2 

2.5)2 

2.5)2 

2.5)2 

2.5) 2 

1.5) 2 
1.5)2 
1.5)2 
1.5)2 


0.1875 

0.1875 

0.1875 

0.1875 

0.1875 

0.1875 

0.1875 

0.1875 

0.1875 


Example 4.4. 
Example 4.5. 


u = 


u = tanh ^ — 30[z — 0.5 — 0.25 sin(27rx) sin(7ry)]^. 
tanh ^ — 3 o|2 ; — tanh ( — 30 [y — 0.5 — 0.25 sin(27rx)]) . 


Adaptive mesh examples (slice and clip cuts) and numerical results are given in Figs. |^to[^ 
As in two dimensions, all three functionals provide good size and shape adaptation, with Q^q 
and Qaii being reasonably small. The best mesh size control Qeq is given by HR functional 


(Figs. 4g, 5g and 6g), although for the considered examples, HR has a slightly worse mesh 


alignment quality Qaii than the others (Figs. 4h, 5h and 6h|. 

A closer look at the example meshes (slice cuts) reveals that, as in 2D, W functional —based 
on variable diffusion— is noticeably more aggressive in moving nodes toward the steep features 
or, alternatively, one can say that the functionals © and ( [I^ based on equidistribution and 
alignment distribute the nodes with the better correspondence with the given M. For coarse 
meshes, all three functionals provide similar results (see convergence plots in Figs. [gI] ); 

however, for fine meshes, sizing of mesh elements obtained by means of W functional is not 
quite as good as for H and HR functionals, as indicated by a larger Qeq. 

Altogether, the linear interpolation error (Figs. and suggests that HR functional 

provides the best mesh, followed by H and W functionals. One may notice from Figs. 
and 1^ that the convergence of the linear interpolation error for W functional slows down near 
= 10^ for Examples 


4.4 


and 4.5, although it seems to improve as the mesh is refined (Fig. |6i|. 


The reason for this behaviour is not clear to us. On the other hand, W functional has the 
simplest form and seems to be more economic to compute than the other two. From tentative 
comparison, mesh generation using W functional uses about one fifth to an half of the CPU 
time used with H or HR functional. Qualitatively, this is not difficult to understand since W 
functional is convex whereas the others are not (although they are polyconvex). 
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(d) Winslow’s (|6| 


10 ^ 10 '‘ 10 ® 10 ® 
(g) Qeq vs. N 


(e) Huang’s (11) 


- Winslow 
Hiiang 
HR 


10 ® 10 '‘ 10 ® 10 ® 
(h) Qah vs. N 


(f) Huang and Russell’s (13) 
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—1— Winslow 
- Hnanff 
HR 
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(i) error vs. N 


Figure 4: Example 4.3 


The top row: slice cuts of the meshes. The middle row: clip cuts of 
The bottom row: Qeq, Qaii, and the L?' norm of the linear interpolation 


the meshes 
error (black line represents 
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(d) Winslow’s ^ 


(e) Huang’s (11) 


(f) Huang and Russell’s (13) 





Figure 5: Example 


4.4 


The top row: slice cuts of the meshes. The middle row: clip cuts of 
the meshes. The bottom row: Qeq, Qaii, and the norm of the linear interpolation 
error (black line represents 
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(a) Winslow’s ^6^ 


(b) Huang’s (111 


(c) Huang and Russell’s (131 






Figure 6 : Example 


4.5, The top row: slice cuts of the meshes. The middle row: clip cuts of 
The bottom row: Qeq, Qaii-, and the L?' norm of the linear interpolation 


the meshes 
error (black line represents 
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5 Conclusions 


Among the three functionals in this study, Huang and Russell’s functional consistently provides 
the best mesh for piecewise linear interpolation in both two and three dimensions. In all examples 
it leads to the best equidistribution quality and the smallest interpolation error. Interestingly, 
while it results in the best mesh alignment quality in two dimensions, the functional gives a 
slightly worse mesh alignment than the other two functionals in three dimensions. 

Meshes obtained by means of Winslow’s functional have the worst mesh equidistribution 
(element size) quality and the largest interpolation error in four out of hve examples, although 
in three dimensions its mesh alignment is quite good and even better than that of the meshes 
obtained using Huang and Russell’s functional. An explanation to this behavior could be the 
fact that this functional does not have an explicit mechanism to control the equidistribution 
property. 

The behavior of Huang’s functional is somewhere in between Winslow’s and Huang and 
Russell’s functionals: both in mesh quality measures and interpolation error. It provides better 
mesh sizing than Winslow’s functional but not quite as good as Huang and Russell’s. On the 
other hand, it provides the best (or very close to the best) mesh alignment in all examples. 

While being able to produce correct and good quality mesh concentration, Winslow’s func¬ 
tional seems to have the tendency to move more points toward the area of interest and is 
slightly less reliable than the other two functionals especially when the mesh is fine. On the 
other hand, it has a very simple form and is more economic to compute than the others. It can 
be a good choice for mesh adaptation at least for coarser meshes, for which all of the three 
functionals produce comparable meshes. 

Finally, it should be pointed out that the numerical experiment we conducted in this work is 
limited and more work is needed to have an extensive and more complete understanding of 
the behavior of the meshing functionals especially in three dimensions. Moreover, the newly 
developed implementation of the variational methods in m has been crucial to the current 
study to perform substantial computations in two and three dimensions. It is our hope that it 
can serve as an efficient tool for use in future studies of mesh adaptation and movement. 
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